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ABSTRACT 

The Orion Nebula Cluster (ONC) is the nearest region of massive star formation and thus a crucial 
testing ground for theoretical models. Of particular interest amongst the ONC's ^ 1000 members 
are: 6*^ Ori C, the most massive binary in the cluster with stars of masses 38 and 9M0 (Kraus et al. 
2009); the Becklin-Neugebauer (EN) object, a 30km s"^ runaway star of ~ SMg (Tan 2004); and 
the Kleinmann-Low (KL) nebula protostar, a highly-obscured, ~ 15 M© object still accreting gas 
while also driving a powerful, apparently "explosive" outflow (Allen & Burton 1993). The unusual 
behavior of BN and KL is much debated: How did BN acquire its high velocity? How is this related 
to massive star formation in the KL nebula? Here we report the results of a systematic survey using 
~ 10^ numerical experiments of gravitational interactions of the O^C and BN stars. We show that 
dynamical ejection of BN from this triple system at its observed velocity leaves behind a binary with 
total energy and eccentricity matching those observed for O^C Five other observed properties of 
6^C are also consistent with it having ejected BN and altogether we estimate there is only a < 10~^ 
probability that 6^G has these properties by chance. We conclude that BN was dynamically ejected 
from the 6*^0 system about 4,500 years ago. BN has then plowed through the KL massive-star- forming 
core within the last 1,000 years causing its recently-enhanced accretion and outflow activity. 
Subject headings: binaries: general - Methods: numerical - Scattering - Stars: individual: Becklin- 
Neugebauer object, 6^G - Stars: kinematics and dynamics 



1. INTRODUCTION 

Massive stars impact many areas of astrophysics. In 
most galactic environments they dominate the radia- 
tive, mechanical and chemical feedback on the interstel- 
lar medium, thus regulating the evolution of galaxies. 
Many low-mass stars form in clusters near massive stars, 
and their protoplanetary disks can be affected by this 
feedback also. There is some evidence that our own so- 
lar system was influenced in this way (e.g., Tachibana 
k. Huss 2003; Adams 2010). Despite this importance, 
there is no consensus on the basic formation mechanism 
of massive stars. Theories range from scaled-up versions 
of low- mass star formation (McKee & Tan 2003), to com- 
petitive Bondi-Hoyle accretion at the center of forming 
star clusters (Bonnell et al. 2001; Wang et al. 2010), to 
stellar collisions (Bonnell et al. 1998). The Orion Nebula 
Cluster (ONC) is the nearest region of massive star for- 
mation and thus a crucial testing ground for theoretical 
models. 

Of particular interest amongst the ONC's ^ 1000 
members in this regard are: 9^ Ori C, the most mas- 
sive binary in the cluster with stars of masses « 38 and 
9Mq (Kraus et al. 2009); the Becklin-Neugebauer (BN) 
object, a 30kms~^ runaway star of « SAf© (Tan 2004); 
and the Kleinmann-Low (KL) nebula protostar, a highly- 
obscured, about 15 M© object still accreting gas while 
also driving a powerful, apparently "explosive" outflow 
(Allen & Burton 1993). The unusual behavior of BN 
and KL is much debated bearing implications towards 
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massive-star formation theories: How did BN acquire its 
high velocity? How is this related to massive star forma- 
tion in the KL nebula? 

BN, like KL, is heavily obscured by dust so its luminos- 
ity of ~ (5 ± 3) X 10'^ Lq mostly emerges in the infrared 
(Gezari et al. 1998). The above luminosity constrains 
BN's mass to be TObn — 9.3 ± 2.0Mq, assuming it is on 
the zero age main sequence (Tan 2004) . For this estimate 
and throughout the paper we have adopted 414 ± 7 pc 
for the distance to the cluster (Menten et al. 2007). As- 
trometry based on mm and radio observations indicate 
that BN is a runaway star (Plambeck et al. 1995; Tan 
2004), with some recent measurements of its motion in 
the ONC frame of ji-qn = 13.2 ± 1.1 masyr^^ towards 
P.A.BN = -27°. 5 ± 4° (Gomez et al. 2008) and /.jbn == 
13.4 ± 1.1 masyr-i towards P.A.bn = -18°.8 ± 4.6° 
(Goddi et al. 2011) (Figure 1). This corresponds to a ve- 
locity V2D,bn = 25.9±2.2kms-i (Gomez et al. 2008). BN 
has an observed radial (LSR) velocity of -|~21± ^ Ikms^^ 
(Scoville et al. 1983), while the ONC mean is -H8.0kms~i 
(based on a mean heliocentric velocity of about 1000 
ONC stars of -^26.1 kms'^ Fiiresz et al. 2008). Includ- 
ing this -|-13± ~ 1 kms~^ radial velocity with respect 
to the ONC mean, the 3D ONC-frame velocity of BN is 
W3D,BN = 29 ± 3 kms""'^. This is much greater than the 
velocity dispersion of ONC stars, variously inferred to be 
(T3D = 2.4 km s^^ based on the proper motions of ^ 50 
bright {V < 12.5) stars within 30' of the ONC center (van 
Altena et al. 1988), ctsd = 3.8 km s""'^ based on proper 
motions of ~ 900 fainter stars within 15' of the ONC cen- 
ter (Jones & Walker 1988), and 0-30 = 5.4kms~^ based 
on radial velocity measurements (potentially affected by 
motion induced by binarity) of 1100 stars within ~60' of 
the ONC center (Fiircsz et al. 2008). Thus there is little 
doubt that BN is a runaway star, which formed and was 
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then accelerated in the ONC. 

Supernova explosion of one member of a binary can 
lead to the other being ejected at high speeds (Zwicky 
1957). The ONC is too young (most stars arc < 3 Myr 
old; Da Rio et al. 2010) for a supernova to have occurred. 
Nor is there any evidence for a recent supernova. Alter- 
natively, runaway stars can be produced via dynamical 
ejection — a gravitational slingshot — from a triple or 
higher multiple system (Poveda et al. 1967; Hut & Bah- 
call 1983; Gies & Bolton 1986), in which the lowest mass 
member tends to be ejected. Indeed, such dynamical 
ejection of stars naturally happen in dense and young 
clusters (e.g., Gualandris et al. 2004; Fujii & Portegies 
Zwart 2011). Similar ejections have also been discussed 
in the context of some other ONC stars (Gualandris et al. 
2004). Thus BN, having formed in the ONC, should have 
been accelerated via dynamical ejection. The predictions 
of this scenario are very specific: somewhere along BN's 
past trajectory should be a massive binary (or higher or- 
der multiple), with two components likely more massive 
than BN, recoiling in the opposite direction, and, as we 
shall see, with specific orbital properties. 

Two scenarios for the dynamical ejection of BN have 
been proposed. (1) Ejection from the 9^C binary (Tan 
2004): In this scenario the BN star is ejected via a strong 
gravitational scattering interaction (e.g.. Hut & Bahcall 
1983; Gvaramadze & Gualandris 2011) and later plows, 
by chance, through the KL star-forming core to drive 
tidally-enhanced accretion and thus outflow activity. If 
so, a model of formation of the KL massive protostar via 
an ordered collapse of a gas core to a central disk(McKee 
& Tan 2003); similar to how low- mass stars are thought 
to form, is still broadly applicable, though subject to the 
tidal perturbation from BN's fly-by. (2) Ejection from 
the KL (source 7) protostar (Bally & Zinnecker 2005; 
Gomez et al. 2008): Here it is proposed that the KL out- 
flow is related to the disintegration of a forming triple 
system, which ejected BN and produced a binary sug- 
gested to be the radio source I (Menten & Reid 1995). 
This binary has recoiled southwards from the original 
formation site and is now hidden, by chance, behind or 
in the dense gas core near the center of the KL nebula. In 
this scenario the gas from the original formation site, like 
the stars, has also been expelled in this event to form the 
outflow, and the core that formed these massive stars has 
been destroyed. If true, this is a very different formation 
process, and would indicate that chaotic gravitational in- 
teractions between multiple protostars followed by com- 
plete ejection of both stars and gas are intrinsic features 
of massive star formation (Bally & Zinnecker 2005; Bally 
et al. 2011), at least in this case in Orion. 

Figure 1 shows a ncar-IR image of the central region 
of the ONC, including BN, the KL protostar (marked 
by radio source T) and the famous Trapezium stars, of 
which 9^C is the brightest. The past trajectory of BN is 
indicated based on its present motion and assuming no 
acceleration. It goes near KL and the Trapezium stars. 
The high obscuration to KL means that there is little 
direct constraint on the properties of the star(s): for ex- 
ample there is no evidence that it is even a binary. In 
contrast, the properties oiO^C have been measured much 
more precisely (Kraus et al. 2009) and so the scenario of 
ejection of BN via a binary-single strong scattering can 
be tested much more rigorously and is the goal of this 



study. 

In §2 we summarize our methods and numerical cal- 
culations. In §3 we present our key results, which show 
that 6^C has orbital properties expected if it ejected BN. 
In §4 we estimate the probability that d^C has not been 
responsible for BN's ejection and has these and other ob- 
served properties simply by chance. We summarize and 
conclude in §5. 

2. METHODS 

We investigate the scenario of ejection of BN from the 
9^C binary by carrying out calculations of gravitational 
scattering between the three stars. We adopt the central 
values of the observed masses of the 9^C binary members 
(Kraus et al. 2009), mgiQ^ = 38.2 ± 3.6 M0 and mgic^ = 
8.8 ± 1.7 Mq. For BN we adopt the central value of the 
mass estimate ttibn = 8.2±2.8 M©, which is based on the 
observed cluster-frame proper motion of 9^C of Hgic = 
2.3±0.2masyr~^ (van Altena et al. 1988), i.e., W2D,eic — 
4.5 ± 0.4 kms~^, and assuming it is due to recoil from 
ejecting BN. Note that this mass is consistent with that 
inferred from the luminosity of BN, discussed above. The 
error range includes an assumed pre-ejection motion of 
the center of mass of the 3 stars along the ejection axis 
of 0.7 masyr"^, i.e., similar to that observed for other 
bright ONC stars (van Altena et al. 1988). 

We carry out a systematic investigation of the three 
possible strong scattering interactions between a binary 
and a single star. Depending on the initial perturber and 
the binary members these interactions can be divided 
into three types. Type 1 - BN star is the perturber: 
Here BN is initially a single star that interacts with 9^C 
members 6*^01 and 9^C2- We denote this initial configu- 
ration as [9^Ci, 9^C2] BN. The square brackets denote a 
bound binary and the star outside the brackets is a sin- 
gle stellar perturber. Type 2 - 9^C2 is the perturber: 
This can be denoted as [9^Ci, BN] 6*1 C2. Type 3 - 9^Ci 
is the perturber: This can be denoted as [9^C2, BN] 
9^Ci. 

Out of all outcomes of our numerical experiments we 
focus on the ones that could lead to ejection of BN. Case 
1 - BN fly-by: This is a Type 1 interaction and the out- 
come of interest is "preservation" where the initial binary 
members remain unchanged and BN flies by after interac- 
tion with the initial binary. We write this interaction as 
[9^Ci,9^C2] BN — > [9^Ci,9^C2] BN. The arrow points 
from the initial to the final configuration. Case 2 - ejec- 
tion of BN from a binary via exchange with 6*^02: 
[6'iCi,BN] 6*102 — > [6*101, 6*102] BN. Case 3 - ejec- 
tion of BN from a binary via exchange with 6*^01: 
[6iiC2,BN] 6iiCi — > [6*101,61102] BN. The "Types" and 
"Oases" are different in the fact that for the Types we 
only take into account the initial conditions and allow all 
outcomes, whereas, the Oases are more restrictive and 
only considers the outcomes where BN is ejected leaving 
behind a bound binary. 

There are 7 parameters that describe the initial con- 
ditions of each binary-single star interaction and since 
the 3-body problem is chaotic we must sample over the 
expected distributions of these parameters: (1) Eccen- 
tricity of the initial binary, Ci. For each case we investi- 
gate two extreme distributions: (A) Circular, e^ = 0, for 
all systems, which may be expected if the binaries have 
recently formed from a gas disk that has damped out 
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Fig. 1. — Near-infrared (J,H,K) image of the central region of the Orion Nebula Cluster (McCaughrean ct al. 2002) with cluster frame 
proper motions of BN (Gomez et al. 2008; Tan 2008) and d^C (van Altcna et al. 1988) indicated by arrows proportional to the size of 
the motion. The coordinates are relative to the position of source I («(J2000)=05 35 14.5141, <5(J2000)=-05 22 30.556) (Gomez ct al. 
2008). Under the simplifying assumption of no acceleration, these motions are traced back with dashed lines (dotted lines indicate la- 
uncertainties) to a common origin about 4,500 years ago, shown by the cross at (-23.72, -44.16). The positions of the other Trapezium 
stars and the Kleinmann-Low massive protostar, source /, are also indicated. 



noncircular motions; (B) Thermal (Heggie & Hut 2003), 
dFb/dei = 2ei, where Fh is the fraction of the binary pop- 
ulation. This is an extreme scenario that would result if 
binaries have had time to thermalize via stellar interac- 
tions with other cluster stars. The actual situation for 
ONC binaries should be between these limits. (2) Semi- 
major axis of the initial binary, a,;. We assume a flat dis- 
tribution (Eggleton et al. 1989) dFb/dlogai = 0.208 from 
0.1 AU (approximately the limit resulting from physi- 
cal contact) to 6300 AU (the hard-soft boundary beyond 
which binaries are expected to be disrupted by inter- 
actions with other cluster stars; Heggie & Hut 2003). 
(3) Initial impact parameter, bi. For each Case and 
each sampling of a^ we investigate the full range of im- 
pact parameters that can lead to scattering events strong 
enough to eject BN with its observed high velocity. This 



is achieved by increasing 6,; from small values until the 
regime where all interactions are weak fly-bys incapable 
of increasing BN's velocity to the observed large value. 
(4) Initial relative velocity at infinity, Vi, in the frame 
of the center of mass of the binary. We assume that 
the stars have velocities following a Maxwellian distribu- 
tion with a dispersion of tT3D = 3 km s~^ (Fiiresz et al. 
2008). We have repeated the numerical experiments with 
(T3D = 2kms^^, finding qualitatively similar results. (5) 
The initial angle of the orbital angular momentum vec- 
tor (Li) of the binary with respect to the velocity vector 
(vi) of the approaching single star, which is assumed to 
be randomly oriented. (6) The initial angle between the 
major axis of the binary orbit and the velocity vector 
of the approaching single star, which is assumed to be 



Chattcrjcc ct al. 



randomly oriented. (7) The initial orbital phase of the 
binary, which is assumed to be random. 

We find cross-sections (E) of the various outcomes of 
the Cases 1, 2, & 3 binary-single interactions numerically 
using the Fewbody software (Fregeau et al. 2004), which 
uses an order 8 Runge-Kutta integrator, by performing 
~ 10^ numerical scattering experiments to sample the 7 
dimensional parameter space that is needed to describe 
the possible interactions. This large number of numerical 
scattering experiments gives us rigorous sampling of all 
properties in the dynamical scattering problem, includ- 
ing the initial semimajor axis of the binary (a^), initial 
eccentricity (e,), binary orbital phases, initial velocity at 
infinity (vi) of the single star, the initial impact param- 
eter of the encounter (bi), the angle between the initial 
major axis relative to Vi and the initial orientation of 
the binary (Li.Vi). For example, we sample bi as fine as 
10~^6oi where bo is the impact parameter at infinity that 
results in a closest approach within 2ai. Starting from 
a small value, hi is sampled with the above-mentioned 
resolution up to at least 6i.,„ax = ^o- Within this in- 
terval smaller intervals of bi, 5bi — lO^^'&o are chosen. 
The impact parameter is chosen from each of these in- 
tervals uniform in the area of the annulus between b^ 
and &,j -|- Sbi. If a particular final outcome of interest 
or "event" is achieved (in particular BN-Velocity or 
BN-True events, defined below), then contribution of that 

event to the total S is simply 6T, = 2TTb^6bi (McMillan 
& Hut 1996; but see more recently Fregeau et al. 2006). 
Afterwards, to ascertain that all energetic encounters are 
sampled, we increase the maximum bi geometrically un- 
til 6i,max = lOO&o- For our case, this large value of &i,max 
corresponds to as large a physical distance as the cluster 
size making sure that all possible energetic encounters 
are captured in the determination of the cross-sections 
E. 

Using the large ensemble of numerical gravitational 
scattering experiments we evaluate the Es for outcomes 
where BN is ejected leaving behind 9^Ci and 9^C2 in 
a bound binary (henceforth, "BN-Ejection" events). A 
subset of the BN-Ejection events where BN is ejected 
with the observed velocity of 29 ± 3 km s~^ are called 
"BN-Velocity" events. We do not put any constraints 
on the binary properties that is left behind for the 
BN-Velocity events. We further calculate the Es of a 
subset of BN-Velocity events where the final binary is 
left with orbital properties similar to those observed of 
0^C, namely, a = 18.13 ± 1.28 AU, and e = 0.592 ± 0.07 
(Kraus et al. 2009, henceforth, "BN-True" events). The 
BN-Ejection events are used to explore the velocity 
distribution of BN if it is ejected via a strong binary- 
single interaction. BN-Velocity events, a subset of the 
BN-Ejection events, show us all possible interactions 
over a range of a^ where BN could have an energy com- 
patible with the observed energy. A further subset, the 
BN-True events, give us stronger constraints and indi- 
cates the range of initial binary properties most likely to 
create the observed ^-'^C binary as well as the runaway 
BN star. 

3. RESULTS 

In this section we present the key results of our numer- 
ical experiments. We start with overall outcomes of all 



our numerical experiments for all Cases and e^ distribu- 
tions and then increasingly focus our attention towards 
the observed BN-6'^C system and compare its various 
properties with those predicted from our simulations. 

Figure 2 shows the branching ratios of all outcomes 
in general from our numerical experiments. A handful 
of interesting aspects are evident in the branching ra- 
tios for the given masses of the 3 stars involved in these 
interactions. 

Disruption (or ionization) of the initial binary happens 
only when the binary is dynamically soft (Hut & Bahcall 
1983), i.e., the value of the binary binding energy is < the 
kinetic energy of the perturber. This is achieved at large 
a?: ^ 300 AU, for Types 1 and 2. For Type 3 ionization 
can happen at relatively smaller a^ '^ 50 AU due to the 
higher mass of the perturber and relatively lower binding 
energy of the initial binary. Nevertheless, even for Type 
3, branching ratio for ionization becomes comparable or 
greater than exchange outcomes only at Ui ^ 10'^ AU. 

For Types 1 and 2, exchange with the primary is very 
unlikely, since here the primary is significantly more mas- 
sive (38.2 Mq) than the secondary (8.8 and 8.2 Mq for 
interactions of Types 1 and 2, respectively). In inter- 
action of Type 3, since the initial binary consists of two 
stars with comparable masses, both exchanges are almost 
equally likely. For Type 3 the fraction of exchange out- 
comes is comparable to the fraction of fly-by events for 
a large range of at taking into account sufficiently strong 
encounters (see §2 for the value of the maximum impact 
parameter). The fraction of fly- by outcomes is of course 
formally infinite since one can always use a sufficiently 
large impact parameter where nothing but a weak fiy-by 
is the outcome. 

Collisional outcomes are comparable with exchanges 
only for sufficiently small a^ values. In interactions of 
Type 1, preservation is the channel that can produce the 
observed BN-6'^C system. For Type 1 for both Ci distri- 
butions collisions become important for fli ^ 1 AU. For 
interactions of Types 2 and 3, exchange of the perturber 
with the primary is the channel that can produce the 
observed BN-0^C system. For Type 2, collisions become 
comparable with the BN-0^C producing channel for Ui < 
a few AU. Whereas, for interactions of Type 3, collisions 
remain comparable to the BN-6'^C producing channel for 
Oi < 1 AU. In interactions of all Types collisions happen 
more often for the Thermal e^ distribution since for the 
Thermal e^ distribution the pericenter distances for the 
stars in binary can be much smaller than that for the 
Circular e^ distribution for any given Oj. We later show 
(§3.3) that BN-True events happen in a given range of a^ 
for a given interaction Type. For all interaction Types 
and Ci distributions collisions have much lower branching 
ratios compared to the branching ratios for the BN-6'^C 
production channels for the ranges of a^ where BN-True 
events can occur. 

In the following sections we increasingly focus on out- 
comes that are similar to the observed BN-6'^C system. 
We first present results for all events where the BN star is 
ejected (BN-Ejection, §3.1). Then we present results for 
all outcomes where the BN star is ejected with a velocity 
within the observed range of 29±3 kms^^ (BN-Velocity, 
§3.2). We then restrict our attention to only a subset of 
the BN-Velocity events where the final binary has prop- 
erties similar to the observed 6^G binary (BN-True, §3.3 



g 
■5 



o 

c 

CO 



-.0 



Runaway star via binary-single interactions 
Type1 Type2 



Types 



"l E™1 — I I iiiiiil — I I iiiiiii I I I r- 



10-' k 

g 

'•t—> 
CO 

^ 10"^ I- 

o 

c 

CO 
^_ 
-Q .r> 

10"* ^ 



10-^ 

10° 



10 ' r 



^ ^0''' r 



10"' 



10-^ r 









^r-- 



.!, iK 



.1 



*'-'" 



.'■■./ 

r 
/\ 

/' 
i 

i " 
jiiiii 



c: 



,'■■ ' -- 



j%1 'r 'v; 




J 



' ! ( I 

I i I 

' I ' I 
I ' I 



^^^\/ 



/:■ 



\ i ' 
/ \ 

i 
i 
i 
■■■■I 



■i I 



"3 i5 
o 

b 






I 

—LUiil l_l-l 



10""' 10° 10'' 10^ 10^10""' 10° 10"' 10^ 10^10""' 10° 10^ 10^ 10^ 



a; (AU) 



aj (AU) 



a; (AU) 



Fig. 2. — Branching ratios for various outcomes in our simulations as a function of a;. Solid (black), dotted (magenta), short-dashed (blue), 
long-dashed (red), and dash-dotted (green) lines denote branching ratios for preservation, coUisional outcome, exchange of the perturber 
with the secondary of the initial binary, exchange of the perturber with the initial primary, and disruption of the binary, respectively. The 
top and bottom panels are for Circular and Thermal e^ distributions, respectively. In both top and bottom panels three panels from left 
to right denote Types 1, 2, and 3, respectively. 



and §3.4). 

3.1. Velocity Distribution of the Ejected BN Star 

Increasing the kinetic energy of BN by about two or- 
ders of magnitude compared to the value expected given 
the ONC's velocity dispersion is at the heart of the prob- 
lem. Hence, we focus on the velocity distribution of BN 
following ejection. In addition, we focus on energy con- 
siderations for the scattering problem, especially the ra- 
tio of the kinetic energy of BN's ejection to the total 
energy of the binary left behind. 

First we explore given the masses of the three stars 
in the interaction, and given that BN is ejected leaving 
the other stars in a binary, how likely it is for BN to 
acquire a velocity significantly higher than the velocity 
dispersion {va = 3km s""'^) in the ONC. We calculate 
the cross-section SBN-Ejection for BN-Ejection events for 
a given a^, for each Cases 1-3, and each e^ distribu- 



tion. The overall cross-section for BN-Ejection events 
for any Oi is calculated by multiplying ^BN-Ejection with 
the probability of finding an initial binary with that a,; 
assuming the semimajor axis distribution for binaries 
is flat in logarithmic intervals within the physical lim- 
its discussed in §2. Thus, the normalized ^BN-Ejection 
is /P(ai)EBN-Ejection(ai,WBN,e)dai, where, P{ai)dai = 
dlog(a,)/log(6310/0.1). 

Figure 3 shows the cumulative distribution of 
SBN-Ejection(i'BN) as a fuuctiou of BN's velocity (vbn) 

calculated using /J^^^2i,„ dweNxd^BN-Ejection/c^UBN- We 
find that binary-single interactions involving the three 
stars in question can eject the BN star with velocities 
that can exceed Vr, by more than two orders of magni- 
tude. However, the cross-sections for such events reduce 
as ubn increases. For both e,; distributions Case 3 shows 
a higher fraction of high-velocity ejection events. This 
is because in Case 3 the perturber is the most massive 
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Fig. 3. — Cumulative distribution of S as a function of the velocity of the ejected star for all cases where BN is ejected. The velocities are 
given in units of the velocity dispersion (vcr = 3kms~^). Top panel is for the Circular e^ distribution. The bottom panel is the same but for 
the Thermal Ci distribution. The solid (black), dashed (red), and dotted (blue) lines in both panels denote Cases 1, 2, and 3, respectively. 
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For each case the cross-section is calculated using f" „ , . ,_,, 

P{ai)dai = 51og(ai)/log(6310/0.1) is used (see text). We use a cut-ofF for the ejection velocity at v 

becomes very large simply due to distant fly- by interactions in Case 1. Note that for all cases dynamical interactions can increase the 

velocity of the ejected star by large factors relative to the velocity dispersion. 



for all events where BN is ejected. For all cases 
2vc, since for v ^ Va S(> v/v^) 



star in the triplet, the one that finally becomes the 6*^01 
star. Hence, the total available energy is higher in Case 
3 events. 

The problem of binary-single scattering can be un- 
derstood by comparing the kinetic energy and the po- 
tential energy of the systems since the outcomes differ 
qualitatively depending on the relative values of these 
quantities (e.g.. Hut & Bahcall 1983). We explore if 
BN is ejected with wbn > 2W(j, then what is the dis- 
tribution of cross-section for the various Cases and e^- 
distributions as a function of the i^ratio and the velocity 
of ejection, wbn- Here, Erratic = Tcjcction/I-Ebinaryl is the 
ratio of the final kinetic energy (Tcjoction) of both the sin- 
gle star and the binary star system (based on the motion 
of its center of mass) to the total energy (gravitational 



energy plus kinetic energy of orbital motion) of the bi- 
nary (E'binary)- Figure 4 shows a 2D distribution of the 
overall SBN-Ejection for any a^ as a function of v^n and 
£^ratio for all Cases 1-3, and all e^ distributions. The 
highest ejection velocities for the BN star happens for 
-Eratio between 0.1 and 1 for all Cases. Note that the 
distributions for Cases 1 and 2 are very similar. This is 
due to the similarity in masses of BN and 6'^C2. How- 
ever, in Case 3 the most massive star is the perturber. In 
addition, the most massive star exchanges into the final 
binary increasing the binding energy of the final binary 
star system. Hence, there is a tail for high i^ratio values 
where BN can still be ejected with a high velocity. Even 
for the more energetic Case 3 events, the maximum ve- 
locity for ejections are achieved within a narrow range of 
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Fig. 4. — 2D distribution of E as a function of energy ratio -Eratio = ^cjcction/l-^'binaryl a-nd velocity (v) of the runaway star BN 
for all events where BN is ejected. The top panels are for the Circular e^ distributions for Cases 1-3 (left to right, respectively). The 
bottom panels are for the same, but for the Thermal ei distribution. The colors denote dE / dE^^i^^ / dv . For each bin the total cross- 
section is calculated using f P{ai)'S{ai, E,-i^tio,s)<io,i for the events where BN is ejected satisfying the ranges of -Eratio and v in that bin. 
P{ai)dai = <51og(ai)/log(6310/0.1) is used (see text). For both e^ distributions CascS shows higher values of -Eratio than Cases 1 and 2. 
This is because in Case 3 the perturber is a more massive star (the one that will finally become 9^Ci ). For all Cases the highest velocity 
increase for the runaway star is between Eratio =0.1 and 1. 



0.1 < -Eratio < 1- 

A high value (~ 10) of the ratio between the veloci- 
ties of the incoming single star (expected to be near Vcr) 
and the ejected single star is needed to create runaway 
stars by definition. This is possible for strong encounters 
involving a binary and a single star where a fraction of 
binding energy of the binary is converted into the kinetic 
energy (T) of the ejected star. If the final outcome is 
again a binary and the runaway ejected star (the binary 
members are not required to remain the initial ones, e.g., 
for Cases 2 and 3) then the kinetic energy of the stars un- 
dergoing dynamical ejection, Tcjoction = (l/2)mBNwiN + 
(l/2)(meici + meic^H^c ^ (6-9 ± 2.7) x 10^6 erg + 
(1.2 ± 0.5) X lO^'s erg = (8.1 ± 2.8) x 10'^^ gj.g (j^gj.g _^ 
indicates the observed values) is expected to be less than 
the magnitude of the total energy of the resulting binary, 
l-Bbinaryl = Gmgic, mgic, /(2a) ^ (16.4±4.9) x 10^6 erg. 
In addition to requiring ii^ratio < 1, one also expects it 
to achieve a value of order unity, i.e. not too much less 
than one. For -Eratio ^ 1, collisional outcomes dominate 
(Figure 2; also see e.g.. Hut & Bahcall 1983; Fregeau 
et al. 2004). On the other hand, for -Eratio ^ 1 outcomes 
with a disruption of the binary dominates creating 3 sin- 
gle stars (Figure 2). The observed value for the BN-0^C 



system is -Eratio — >■ 0.49 ±0.22, consistent with it being a 
result of a binary-single interaction. 

3.2. Kinetic Energy of Ejection and Eccentricity of 
Binary 

We now focus on the subset of BN-Ejection events 
(BN-Velocity) where BN is ejected with the observed 
velocity of 29 ± 3 (§2) via any one of the Cases 1-3 (a 
small strip in the vertical axis in Figure 4). The 6^C 
orbit is eccentric (e « 0.6). Indeed, strong encounters 
are expected to leave behind binaries with generally high 
eccentricities. We now want to see the 2D distribution of 
overall cross-section for all BN-Velocity events for any 
Oi as a function of -Eratio and the final eccentricity. 

In Figure 5 we plot the 2D distributions of S for the 
BN-Velocity events in the -Eratio versus final e plane for 
Cases 1, 2, & 3 for both the Circular and Thermal e^ dis- 
tributions. The observed values of the BN-0^C system 
are also shown. These overlap within the la contours 
for all cases. The energy of the 9^C binary is just what 
we would expect if it had ejected BN at the observed 
velocity. Its high eccentricity, ~ 0.6, is also naturally ex- 
plained by the recent ejection of BN since during ejection 
the potential of the system is changing rapidly. 

Note that if the d^C binary was unrelated to BN, then 
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Fig. 5. — 2D distribution of S as a function of energy ratio E,.; 
the BN-Velocity events for all cases and e^ distributions. The top panels are for the Circular e; distributions for Cases 1-3 (left to right, 
respectively). The bottom panels are for the same, but for the Thermal e; distribution. The colors denote dYj j dE-^^XAo I de. For each bin the 
total cross-section is calculated using / P{aiyi^(ai, E^atioj 6)'^'ii for the BN-Velocity events satisfying the ranges of i?ratio a^nd £ in that bin. 
P{ai)dai = (51og(ai)/log(6310/0.1) is used (see text). The point and errorbars show the observed BN-S^C system. The larger horizontal 
errorbars (red) denote the -Eratio errors including contribution from the mass errors of the stars, whereas the shorter (white) errorbars 
denote the same if the central mass values are chosen and no contribution from mass errors are included (consistent with the numerical 
experiments). For all cases and all e^-distributions the observed system properties lie within the la contours of the E-distribution, with 
Case 3 Circular e^ being somewhat more favored (see text). 



£^ratio could liave values in a large range spanning many 
orders of magnitude. For example, the range in a for the 
binary is determined by contact (^ a few stellar radii) 
to the hard-soft boundary in ONC {^ 6000 AU for the 
velocity dispersion in ONC). Thus i^binary for 9^C if un- 
related to BN's ejection, could be expected to have values 
anywhere between ^ 10^'* and lO"'^ ergs. It is thus very 
interesting to find the observed BN-6'"'^C system with en- 
ergies so close to the expected energies if they have had a 
binary-single scattering encounter in the past. The like- 
lihood of this occuring simply by chance rather than by 
being caused by interaction with BN is explored in §4. 

3.3. Cross-Sections o/ BN-Velocity and BN-True 
Events: 

We evaluate the cross-sections for outcomes where BN 
is ejected with the observed velocity of 29 ± 3km s^^ 
("BN-Velocity" events). We further calculate the cross- 
section of a subset of BN-Velocity events where the fi- 
nal binary is left with orbital properties similar to those 
observed of O^C, namely, a = 18.13 ± 1.28 AU, and 
e = 0.592 ± 0.07 (Kraus et al. 2009; "BN-True" events). 

For Cases 1, 2, 3 with Circular initial e; distribu- 
tion, SBN-Vciocity = (2.7,0.94,2.9) X 103 AU^ while 



for Thermal initial e^ distribution, SBN-Vciocity = 
(2.9,1.4,2.5) X 10^ AU^ For Cases 1, 2, 3 with Cir- 
cular ei-distribution, Sbn-Ttuo = 56, 24, 105 AU , while 
for thermal e^ distribution, Sbn-Tiuc = 58,34, 84 AU'^. 
The three cases have similar cross-sections, with Case 
3, i.e. [6'iC2,BN] 6*101 — > [6liCi, 6*102] BN, somewhat 
more preferred. 

The cross-sections for BN-Velocity and BN-True 
events as a function of bi, are shown in Figure 6. The 
cross-sections for a given bi and for all explored a^ are 
normalized using the probability, P{ai), of finding a bi- 
nary with semimajor axis a^, assuming a semimajor axis 
distribution flat in log intervals. For small values of 
bi < 300 AU, the cross-sections grow geometrically as 
hj and then decline. This decline, seen in all cases, is 
simply due to the fact that the interactions are happen- 
ing at larger and larger impact parameters and at some 
point no interactions are expected to be strong enough 
to increase the energy of the ejected star to the observed 
value of BN. The areas under the histograms are the 
total cross-sections for the BN-Velocity and BN-True 
events. Figure 6 also includes a table summarizing the 
total cross-sections. 
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Fig. 6. — Cross-section, S, vs initial impact parameter, 6^, of interaction. The top and bottom panels show results for the Circular 
and Thermal a distributions, respectively. Thin and thick lines represent BN-Velocity and BN-True events, respectively. Black (solid), 
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denote the total S values for BN-Veloclty and BN-True events for all cases and ei-distributions. All three cases can contribute towards 
producing the observed BN-9^C system, however. Case 3 is preferred by a factor of ~ 2. 



The cross-sections for BN-Velocity and BN-True 
events as a function of a^ are shown in Figure 7. Wc 
find that BN-Velocity events happen via binary-single 
scattering encounters for a large range of a^, but BN-True 
events place much tighter constraints on a^. For example, 
Case 3, which is the most favored, requires the initial bi- 
nary [0^C2, BN], i.e. two approximately equal-mass stars 
of ~ 9Mq, to have originally had a^ ~ 8 ± 2 AU. Note 
that depending on the case, different ranges of a^ con- 
tribute towards creating systems similar to the observed 
BN-6'^C system. Moreover, note that the ai ranges where 
the BN-True events occur via the three Cases 1 ~ 3, ejec- 
tion of BN is the most likely outcome among all other 
possible outcomes (except of course weak fly-bys; Fig- 
ure 2). This gives us further confidence in the scenario 
that the BN-0^C system has been created via a strong 
binary-single interaction involving these three stars. 

Throughout this study we have used the central values 



of the estimated masses for the three stars. Dynamically 
there should be no qualitative difference in the outcomes 
if the masses are changed within the mass errors. How- 
ever, note that the masses of 9^C2 and BN are compa- 
rable. In fact the error ranges actually overlap. Due 
to the comparable masses, dynamically there is only a 
small difference between ejection of BN and ejection of 
0^C2, as named here. This is reflected in our results to 
some degree. Cases 1 and 2 contribute towards creating 
the BN-0^C system over very similar a^ ranges. Their 
contributions are also comparable. The small differences 
in the ai range and ^BN-True come from the small differ- 
ence between 9^C2 and BN's assumed masses and also 
to some extent the scenario, a fly-by being more likely 
compared to an exchange if all else is kept unchanged. In 
fact, similar results will be recovered if BN and 9^C2 are 
interchanged among themselves. However, in that case, 
the definitions of Cases 1 and 2 will also be interchanged. 
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Fig. 7. — Same as Figure 3, but for S as a function of the initial seniimajor axis, a^, of the binary. Wc find that all three cases and 
ei-distributions can contribute to the production of the observed system over different Oi-rangcs, although Case 3 has a relatively higher 
integrated S for both ei-distributions and selects from a quite distinct range of a^ ~ 8 ± 2 AU. 



3.4. Orientation of the Orbital Plane of 9^G Binary 
Relative to the Direction of BN's Motion 

One additional variable in the dynamical ejection prob- 
lem is the angle a between the angular momentum vec- 
tor, Zflic, of the O^C binary and the 3-D velocity vector 
of BN, wbn- Figure 8 shows the distribution of SsN-True 
for all BN-True events for any Oi as a function of a for 
all Cases 1-3 and all Ci distributions. The distributions 



of a are quite broad for Cases 1 and 2. In comparison, 
for Case 3 there is a strong peak near a = 90°. 

To calculate the observed value of a for comparison 
with the predictions of our numerical results we adopt 
different measured values in existing literature. Orien- 
tation of LgiQ is calculated using data given in Kraus 
et al. (2009). The LSR velocity of BN is obtained from 
Scoville et al. (1983). There are two independent proper- 
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Orientation of the 9^C binary is obtained from Kraus et al. (2009). The predicted distributions from our simulations are consistent with 
the estimated values of a for the observed BN-0^C system. 



motion measurements for BN. Adopting the values given 
in Gomez et al. (2008) we find a = 95°. Adopting the 
values given in Goddi et al. (2011) we find a slightly dif- 
ferent value of a = 87°. If the direction of BN's motion 
in the sky-plane is obtained by simply joining the ex- 
pected position of the binary-single encounter and BN's 
current position, then a = 99°. (this may be a more 
accurate value, since we expect BN to have suffered a 
recent change in its proper motion vector via interaction 
with source /, see below) Note all of the above values of 
a for the observed BN-6'^C system are consistent with 
the predicted distribution of a from our numerical ex- 
periments, and give some support for the ejection having 
resulted via Case 3. 

4. TESTS FOR THE EJECTION SCENARIO OF BN FROM 
e^C AND PROBABILITY OF CHANCE AGREEMENT 



The system which ejected BN must be located along 
BN's past trajectory and have a total mass > 2tobn- 
These conditions are potentially satisfied for d^C, the 
3 other Trapezium stars O'^A, 6^B, d'^B, another ONC 
member 9'^ A, and probably for source / (assuming it is 
the main source of luminosity in the KL nebula). Indeed, 
a number of authors have argued BN was launched from 
source / (Bally & Zinnecker 2005; Gomez et al. 2008). 
However, as we now discuss, there are 6 independent 
properties of 9^C that have the values expected if it were 
the binary left behind after ejecting BN (7 if we assume 
ejection via Case 3 and include the angle a between the 
angular momentum vector, LgiQ, of the 0^C binary and 
the 3D velocity vector of BN, wbn)- To consider the like- 
lihood that all of these properties of the BN-6'^ C system 
are as observed by chance^ we take that to be our null 
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hypothesis. Given that BN has the runaway velocity, for 
each of these properties we assign a probabiUty that 9^C 
has its values by chance to finally calculate the compos- 
ite probability of chance agreement of fl^C's properties 
with those expected from a binary-single scattering sce- 
nario. We discuss these properties and our estimates of 
the chance- agreement probabilities below. These proba- 
bilities are summarized in Table 1 along with the values 
predicted by the binary-single ejection scenario, and the 
observed values of the BN-6'^C properties. 

(1) ONC- Frame Proper Motion in Declination 
(tis,ONc{6^C)): If 9^C ejected BN, then, in the frame 
of the center of mass of the pre-ejection triple, the pre- 
dicted value of ij,s,t{&^G) = — (mBN/"T-gic)M<5,T(BN) — )- 
-([9.3 ± 2.0Mq]/[47 ± 4Mo])11.7 ± 1.3 masyr'^ -^ 
—2.3 ± 0.6 mas yr~^. Here we have used the luminosity- 
based mass estimate for BN (Tan 2004) and the 
proper motion measurements of Gomez et al. (2008) 
for BN, including a 0.70 mas/yr uncertainty of 
the motion in declination of the pre-ejection triple 
with respect to the ONC frame. The predicted 
value of the ONC-frame motion of 9^C is then 
M(5,ONc(^^C) [predicted] = —2.3 ± 0.9 masyr^^, 
with the error increasing again because of the 
uncertain motion of the pre-ejection triple^. 
The observed value (van Altena et al. 1988) is 
M<5,ONc(^^C') [observed] = — 1.8 ± 0.2 masyr"^. Given 
the observed (van Altena et al. 1988) ID proper motion 
dispersion of the bright ONC stars of 0.7 mas/yr, the 
probability for 6^C to be in the predicted range is 0.023. 
Indeed, van Altena et al. (1988) already noted that 9^C 
has an abnormally large proper motion. 

(2) ONC-Frame Proper Motion in Right Ascension 
(^^a,ONc{S^C)cos6): Similarly, in the frame of the cen- 
ter of mass of the pre-ejection triple: /iQ^T(^^C)cos(5 — 
-(mBN/TOeic)Aia,T(BN)cosJ -^ -([9.3 ±'2.0Mo]/[47 ± 
4M0])(-6.1± 1.2) masyr-i -^ +1.21 ± 0.36 mas yr^^. 
The predicted value of the ONC-frame motion of 
9^C is then /ia^oNc(^^0)cos(5[predictcd] = -1-1.2 ± 
0.8mas yr~^. The observed value (van Altena et al. 1988) 
is /iQ,eicCOS(5[observed] = -1-1.4 ± 0.2 mas yr~^ and the 
probability that 9^C is in the predicted range by chance 
is 0.27. 

(3) ONC-Frame Radial Velocity Cwr,ONc(^^C);.- 
Similarly, the radial recoil in the frame of the 
pre-ejection triple should satisfy: WrT(^^C) = 
-(mBN/m,ic)«r,T(BN) ^ -([9.3 ± 2.6Mo]/[47 ± 
4Mo])(-Hl3 ± 1.8 km s-i) -^ -2.57 ± 0.69 km s^^. 
The predicted value of the ONC-frame motion of 9^C 
is then Wr,ONc(^^C) [predicted] = —2.6 ± 1.6 km s~^, 
with the error range dominated by the assumption 
that the pre-ejection triple had a motion similar to the 
other bright ONC stars (van Altena et al. 1988) with 
(TiD = 1.4 km s~^. 9^C has an observed heliocentric 
velocity (Kraus et al. 2009) of +23.6 km s"! i.e. an LSR 
velocity of 5.5 km s~^, i.e. an ONC frame velocity of 
Wr,ONc(^"'^C) [observed] ~ — 2.5kms~^. For a Gaussian 
distribution with ctid = 1.4kms~^, i.e. based on the 
proper motion dispersion of bright stars (van Altena 

^ Note that here and for the other 9^C properties we have 
adopted la errors when possible, but not all physical properties 
have well-defined uncertainties: e.g., the model dependent niBN 
estimate given its observed luminosity. 



et al. 1988) the probability of being in the predicted 
velocity range by chance is 0.24. 

(4) Mass of Secondary (mgiQ^): Given a 9^C primary 
mass of 38.2 M0, what is the probability of having a 
secondary star with mass > ttibn? We estimate this 
probability using the low value (tobn = 7.3 Mq) of the 
luminosity- based mass estimate for BN of 9.3 ± 2.0 Mq. 
If the secondary star is drawn from a Salpeter power- 
law mass function, dF/drntf ex ttt.^^'^^, where F is the 
fraction of the stellar population, with maximum mass 
equal to the primary mass and lower mass limit equal to 
1.0 M0 (a relatively top-heavy IMF, with average mass 
of 2.8 Mq, compared to the global ONC IMF, which has 
a broad peak around 0.5 M©, e.g., Muench et al. 2002), 
then the probability of having a secondary with mass 
> 7.3 Mq is 0.061. Hillenbrand & Hartmann (1998) find 
evidence for mass segregation in the IMF [of primary 
stars) in the center of the ONC, with average stellar mass 
reaching peak values of ~ 1.3 — 2 Af© in the vicinity of 
9^C (excluding 9^C from the average). If we raise the 
lower limit of the above Salpeter IMF to 2 M© (i.e. an 
average mass of 5.1M©), this raises the probability of ob- 
taining a sufficiently massive secondary to 0.16 and we 
adopt this number as a conservative estimate. We note 
that none of the other Trapezium stars has a secondary 
mass that satisfies this condition (Tan 2004). 

(5) Ratio of Ejection Kinetic Energy to Binary Total 
Energy (E-cs.tio(BN-9^G)): From our numerical experi- 
ments we find that in order for 9^C to have ejected BN 
at its observed velocity, 0.23 < E^ratio [predicted] < 0.72. 
This range contains about 70% of the BN-Velocity 
events. Recall, -Eiatio [observed] — 0.49 ± 0.22. Given 
the primary and secondary masses of ^"'^C, what is the 
probability the total energy of the binary, i^binary, falls 
in the predicted range of 1.4 to 4.3 x Tojoction, i-e. 
(1.1 — 3.5) X 10"*^ erg, simply by chance? This corre- 
sponds to a range of semi-major axes of 8.5 to 27 AU. If 
the distribution of a follows dFb/d\oga = constant (Heg- 
gie & Hut 2003) from 0.1 to 6300 AU (see §2), then 
the probability that the 9^C binary falls in this range 
is 0.10. If we assume ejection occurred via Case 3, then 
0.44 < E'ratio [predicted] < 0.83, corresponding to a range 
of semi-major axes of 16 to 30 AU and a probability of 
chance agreement of 0.057. The upper limit of allowed a 
may be smaller than 6300 AU for conditions in the central 
regions of the ONC. For a stellar density of ^ 10** pc~^, 
the average separation is 6000 AU. Most of these stars 
will have masses lower than 6'^C20r BN. Hence it is not 
likely that these lower mass stars will disrupt the rela- 
tively more massive BN or 9^C2 stars from being in a 
binary with 9^Ci. If we adopt an upper limit smaller 
by a factor of 2, i.e. 3000 AU, then the probabilities of 
the 9^C binary falling in the above expected ranges by 
chance rises by just 7%. 

(6) Eccentricity (e(9^C)): Our numerical experiments 
(see Figure 2) show that a very broad range of eccen- 
tricities is expected for the 6*^0 binary if it has ejected 
BN at the observed velocity (the average value for all 
the BN-Velocity outcomes is e = 0.60 with a \a range 
from 0.34 to 0.86; for Case 3 this range is from 0.34 to 
0.82). The observed value of e = 0.592 ± 0.07 is con- 
sistent with these expectations, especially being close to 
the peak of the distribution resulting from Case 3 (Figure 
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Fig. 9. — Cumulative histogram weighted by respective cross-sections for the closest approach of either BN or 0^02 to Q^C\ {Rmin eiCi ) 
during the binary-single interactions for all cases and ei-distributions for BN-True events. For all cases a significant fraction of interactions 
results in R^i„ fliCi ~ 5AU. Especially, for Case 3 and Circular e^-distribution, probability for R^i„ eiCi ~ 5AU is about 60%. 



5). However, to assess the probability of this agreement 
by chance we need to know the eccentricity distribution 
of ONC binaries, especially for massive stars. Unfortu- 
nately there are few observational constraints on this ec- 
centricity distribution. If the binaries have existed long 
enough to suffer many interactions, then a thermal dis- 
tribution, dFb/de = 2e, is expected, which is weighted 
towards high eccentricities. A binary drawn from a ther- 
mal distribution has a 0.62 chance to be in the range 
0.34 < e < 0.86 (0.56 to be in the range 0.34 < e < 0.82 
for Case 3). Of course, if the actual distribution of e of 
ONC binaries is close to circular (e ~ 0), then the prob- 
ability of chance agreement for the eccentricity oi 9^C 
with the value expected from BN ejection would be very 
small. To be conservative, we adopt the probability of 
0.62 implied by a thermal distribution of eccentricities. 

(7) Angle between the direction of BN's motion and 
the angular momentum of O^C binary (a): From our 
numerical experiments we find that for BN-Velocity 
events the angle a between Lqiq and wen should satisfy 
54° < a < 126°. This range contains about 70% of aU 
BN-True events. For Case 3, the range is 58° < a < 122°. 
The observed value of a for the BN-0^C system is ~ 90° 
and thus contained within this range. If 0^C and BN were 
unrelated, then the distribution of a should be 0.5 sin a 
over the range 0° to 180°. Hence the probability for 
chance occurrence for a to be within the above range is 



0.59 (0.53 for Case 3). 

Combining the above individual probabilities and as- 
suming, reasonably, that these properties are mutually 
independent, we find that the total probability of chance 
agreement of all 7 properties is small, e ^ 0.023 x 0.27 x 
0.24 X 0.16 X 0.10 X 0.62 x 0.59 = 8.7 x lO^^. if we 
assume ejection happened via Case 3, which affects the 
last three probabilities, we obtain a total probability of 
chance agreement ofe = 4. Ox 10~^. The excellent agree- 
ment between the observed BN-^-'^C-system properties 
and the predicted final properties from our numerical 
simulations (§3) together with the low chance agreement 
probability (e) of all independent observed properties of 
the system strongly supports our proposed scenario that 
the observed system resulted from a strong binary-single 
ejection event involving 6*^01, 9^G2, and BN with a high 
probability of about 1 - e « 0.99999. 

From the present day velocities and projected distance 
between the 6*^0 and BN this binary-single encounter 
must have happened about 4500 years ago at the loca- 
tion shown in Fig. 1. This scenario also explains the 
anomalously large proper motion of 9^C (van Altena 
et al. 1988), which will cause it to leave the central region 
of the cluster within ^ 10^ years. 

Based on a multi-frequency radial velocity analysis, 
Lehmann et al. (2010) have suggested 6*^0 may actu- 
ally harbor an additional star of 1.0 ± 0.16 M© in a close 
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(a = 0.98 AU; P = 61.5 day), eccentric (e = 0.49) or- 
bit around 9^Ci. If during the proposed scattering event 
that ejected BN, 9^C2 or BN came close to this inner 
region then one would expect likely ejection of the so- 
lar mass star. To examine the likelihood of such close 
interactions, in Figure 9 we show the cumulative his- 
togram (weighted by cross-section) of closest approaches 
of either 6*^02 or BN to ^^Ci for ah BN-True events. 
For example for Case 3 with Circular e^, ^ 60% of the 
events happen with a closest approach that is > 5 AU. 
At these distances we would expect the solar mass star 
to remain relatively undisturbed in its orbit. Case 1 and 
especially Case 2 involve somewhat closer approaches, 
although both of them still have at least a 20% contri- 
bution to events with closest approach > 5AU. Future 
confirmation of the reality of the third star in the 9^C 
system may help us place further constraints on how BN 
was ejected from 6^C. 

5. DISCUSSION AND SUMMARY 

We present the following argument in favor of our pro- 
posed scenario that BN was ejected by a binary-single 
interaction involving 6'^Ci, 9^C2, and BN in the past. 
BN is a runaway star (Plambeck et al. 1995; Tan 2004). 
It had to be launched by an interaction with a multi- 
ple system that has a primary mass greater than BN's 
mass along its past trajectory. The most massive binary 
in the ONC, 9^C, is a system satisfying this condition, 
but there are a few other candidates including, poten- 
tially, the massive protostar source I. To test whether 9^C 
binary ejected BN we consider 7 additional properties, 
namely recoil in 3 directions, sufficiently massive sec- 
ondary, orbital binding energy, orbital eccentricity and 
angle between the orbital angular momentum vector and 
the direction of BN's velocity. Aided in part by a large 
and well-sampled suite of numerical simulations we show 
that all of these 7 observed properties of the BN-6'^C sys- 
tem agree well with the properties predicted if BN was 
ejected by 9^C. There are two and only two possibilities: 
1) 9^C has all these properties by chance; 2) 9^C has 
acquired these properties naturally as a result of BN's 
ejection. We estimate the probability of chance agree- 
ment for all of the above properties to be low e < 10^^. 
Hence, we conclude with about (1 — e) w (1 — 10"^) 
probability that today's BN-0^C system was created via 
a binary-single interaction involving ^"'^Ci, 9^C2, and BN. 

A summary of our numerical calculations is as fol- 
lows. We performed ^ 10^ numerical simulations that 
allowed us to properly sample the multidimensional pa- 
rameter space effectively (§2). We used two different 
Ci distributions as limiting cases, and 200 different at 
values from the full range of possible a^ in ONC from 
physical considerations for all three cases that can pro- 
duce the observed BN-0^C system (§2). We found that 
the predicted energies of the BN-0^C system agree well 
with the predictions from the scattering scenario (Fig- 
ure 5) for both assumed e,; distributions and all three 
cases. Furthermore, the predicted distribution of the an- 
gle a between BN's velocity vector ubn and fl^C's an- 
gular momentum vector LgiQ was consistent with the 
observed value of a (Figure 8). Further calculations of 
cross-sections for the BN-Velocity and BN-True events 
constrained a range of initial binary semimajor axes for 
each Case (1, 2, 3) that would produce the observed BN- 



9^C system followed by a strong scattering (Figure 7). 
For these ranges of a^ ejection of BN in general is the 
dominant outcome (apart from weak fly-bys which has 
a formally infinite cross-section; Figure 2). Our results 
indicate that all Cases can contribute to the produc- 
tion of the observed BN-6'^C system. However, an in- 
teraction between an initial binary with members 9^C2 
and BN, and a single star 9^Ci, our Case 3 (denoted 
by [6*102, BN] ^^Ci -^ [6*101,61102] BN), is favored by 
about a factor of 2 (Figure 6). 

Ejection of BN from 9^C has several important impli- 
cations for our understanding of massive star formation 
in the KL nebula, where a core of gas appears to be col- 
lapsing to form at least one massive star, thought to be 
detected in the radio as source / (Figure 1). 

If BN was ejected from 9^C, then its passage near 
source / within ^ 0.5", i.e. a projected separation of 
~ 200 AU, i.e. an expected physical separation (Gomez 
et al. 2008; Goddi et al. 2011) of ~ 300 AU in the KL neb- 
ula is coincidental. Our estimated ejection point (Figure 
1) is 50.1" from source Fs present location, corresponding 
to 20,800 AU. The ONC-frame radial velocity is about 
half of the plane of sky velocity, implying BN also trav- 
elled 10,400 AU in the radial direction to reach source 
/, for a total distance of 23,200 AU. The probability to 
approach within 300 AU of source /, ignoring gravita- 
tional focussing, is thus 7r(300/23, 200)2/47r = 4 x lO-^. 
Gravitational focussing by ~ 20 M© of total mass in and 
around source / (McKee & Tan 2003; Wright et al. 1992) 
boosts the cross-section by ^14%, so the probability of 
approach is ^ 5 x 10~^. Thus the interaction of BN with 
KL is an improbable event. Of course, the chance of 
interaction of BN with any existing protostar, not neces- 
sarily source /, in the ONC is larger simply by a factor 
equal to the number of protostars in this volume around 
9^C. From X-ray observations (Grosso et al. 2005) there 
appear to be at least ^10 such objects even in just the 
local vicinity of KL, so the total probability of an interac- 
tion between BN and a protostar can be at least an order 
of magnitude higher, but still leaving it as being highly 
unlikely. Wc note that the probability that 6*10 is mas- 
querading as the system that ejected BN is even smaller 
than the probability of BN interacting with source /. Fur- 
thermore, the evidence for ^^C to have ejected BN is 
based on 7 independent lines of evidence and so survived 
several tests for falsification (§4). 

A close passage of BN with source / will have deflected 
BN's motion by an angle 

7.5°(m7/2OM0)(b/3OOAU)-i(vBN/3Okms-i)-2 to- 

wards source /, where b is the initial impact parameter 
and wbn is the velocity of BN relative to source /. 
Wc expect that a reasonably accurate estimate of the 
original ONC-frame position angle (P.A.) of BN's proper 
motion can be derived by considering the angle from 
the estimated position of the dynamical ejection from 
9^C (the cross in Figure 1) and BN's current position, 
which is -29°. 8. The current observed P.A. of BN's 
ONC-frame proper motion is variously estimated to be 
-27°.5±4°(G6mez et al. 2008) and -18°.8±4.6°(Goddi 
et al. 2011). The average of these is -23°. 1 ± 3°, 
suggesting a projected deflection of 6°.7± ^ 3° towards 
source /. The total true deflection may be expected to 
be ~ ^/2 larger, i.e. ~ 9°.5 ± 3°, which is consistent 
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TABLE 1 
Current properties of O^C that are required if it ejected BN 



Property oi e^C 


Predicted Value 


Observed Value^ 


Probability of 
Chance Agreement 


Proper Motion in Dec.*^ if-S.OKc) 


-2.3 ±0.9 mas yr~' 


-1.8±0.2masyr"i 


0.023 


Proper Motion in R.A.'' (/^q,onccos5) 


+1.2±0.8masyr~' 


+1.4±0.2masyr~l 


0.27 


Radial Velocity'' (iJr.ONc) 


-2.6±1.6kms-i 


-2.5kms-i 


0.24 


Mass of Secondary [mgiQ^) 


7.3- 38.2 Mq 


8.8±1.7M0 


0.16<= 


Eject. KE to Binary Total E (Eratio) 


0.23 - 0.72 


0.49 ±0.22 


0.10 [0.057]<J 


Eccentricity (e) 


0.34 - 0.86 


0.592 ± 0.07 


0.62°[0.56]'* 


Angle between LgiQ & {)bn (a) 


54° - 126° 


~90° 


0.59[0.53]<i 


Combination of 7 Independent Properties 






8.7 X 10-''[4.0 X 10-'']'^ 



''References: 9^C proper motions from van Altena ct al. (1988); Other properties from Kraus ct al. (2009). 

''ONC- frame 

'^Assumes secondary is drawn from a Salpctcr initial mass function (IMF) with lower mass limit of 2 Mq. A lower 
limit of IMq (still a top-heavy IMF) would yield a probability of 0.061. 

•^Number in square brackets assumes ejection via Case 3. 

"Assumes binary eccentricity is drawn from a thermal distribution, which is the most eccentric distribution that 
can be expected (requiring the cluster stars to have had a long enough time to interact). Thus this probability 
should be regarded as a conservative upper limit. 



with our previous estimate. 

Close passage of BN near the accretion disk of source 
/ about 500 years ago would induce tidal perturbations 
that enhance the effective viscosity of the disk and thus 
its accretion rate (Tan 2004). This can explain the en- 
hanced, apparently explosive, outflow (Allen & Burton 
1993), which has observed timescales of ^ 500 — 1000 yr, 
if the inner region of the disk with an orbital time 
< 500 yr has been significantly perturbed. This corre- 
sponds to disk radii of r < 18O(m//2OM0)i/3 AU, which 
is consistent with the estimate of closest approach based 
on deflection of BN's trajectory. 

The estimated current accretion rate to source / is 
~ 3 X 10^'^ M0yr^^ to a 10 M0 protostar (and somewhat 
higher rates if the protostar is more massive) (Testi et al. 
2010). This is about a factor of 10 higher than expected 
given the properties of the gas core (McKee & Tan 2002, 
2003). If constant over the last 1000 yr, this would im- 
ply a total accreted mass of about 3 Mq , which would be 
a significant fraction of the original accretion disk-mass 
around a 10 — 20 M0 protostar, since the disk mass is 
expected to be limited to ~ 30% — 100% of the stellar 
mass by gravitational torques (Kratter et al. 2010). The 
mass launched by a magneto-hydrodynamic outfiow dur- 
ing this time is expected to be « 30% of this amount 
(Najita & Shu 1994), i.e. about 1 Mq. This is consistent 
with the mass estimated (Chernoff et al. 1982) to be in 
the inner, "explosive" part of the outflow of about 3 M0. 

Our results suggest that the formation of a massive 
star in the KL nebula, i.e. source /, has been affected 
by an external perturbation of a runaway B star, BN, 
ejected from a different region of the cluster, i.e. by d^C 
For BN to be launched so close to a forming massive pro- 
tostar does appear to be an intrinsically unlikely event, 
but multiple pieces of independent evidence strongly sup- 
port this scenario. It also explains source Fs anomalously 



high accretion rate and the unusual, apparently "explo- 
sive" nature of the recent outflow from this source. It is 
possible that a significant fraction, up to ^ 1/3, of the 
accreted mass has been induced by this perturbation. In 
other respects, the Core Accretion model (McKee & Tan 
2002, 2003) of massive star formation provides a reason- 
able description of the system, e.g., the presence of a 
massive core around the protostar and two wide-angle 
outflow cavities from which near-IR light emerges (Testi 
et al. 2010). The example of Orion BN-KL suggests that, 
occasionally, in crowded regions near the center of star 
clusters, the star formation model needs to be modifled 
to account for tidal perturbations from external, passing 
stars. 

Improved observational constraints on the properties 
of 9^C especially its binary properties and ONC-frame 
proper motion, and the mass and current proper mo- 
tion of BN will place even more stringent tests on the 
proposed ejection scenario and also help improve the dy- 
namical mass constraints on source / from the deflection 
of BN. 
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